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ABSTRACT 

We present theoretical models of the time-dependent thermal and chemical structure of molecular gas sud- 
denly exposed to far-ultraviolet (FUV) (6 eV < hv < 13.6 eV) radiation fields and the consequent time- 
dependent infrared emission of the gas. We focus on the response of molecular hydrogen for cloud densities 
ranging from n = 10 3 to 10 6 cm 3 and FUV fluxes G 0 = 10 3 -10 6 times the local FUV interstellar flux. For 
G 0 /n > 1(T 2 cm 3 , the emergent H 2 vibrational line intensities are initially larger than the final equilibrium 
values. The H 2 lines are excited by FUV fluorescence and by collisional excitation in warm gas. Most of the 
H 2 intensity is generated at a characteristic hydrogen column density of /V — 10 21 cm 2 , which corresponds 
to an FUV optical depth of unity caused by dust opacity. The time dependence of the H 2 intensities arises 
because the initial abundances of H 2 at these depths is much higher than the equilibrium values, so that H 2 
initially competes more effectively with dust in absorbing FUV photons. Considerable column densities of 
warm ( T - 1000) K H 2 gas can be produced by the FUV pumping of H 2 vibrational levels followed by col- 
lisional de-excitation, which transfers the energy to heat. In dense (n ^ 10 5 cm 3 ) gas exposed to high (G 0 > 
10 4 ) fluxes, this warm gas produces a 2-1 S(l)/l-0 S( 1) H 2 line ratio of -0.1, which mimics the ratio found in 
shocked gas. In lower density regions, the FUV pumping produces a pure-fluorescent ratio of -0.5. We also 
present calculations of the time dependence of the atomic hydrogen column densities and of the intensities of 
O i 6300 A, S n 6730 A, Fe n 1.64 gm, and rotational OH and H 2 0 emission. Potential applications include 
star-forming regions, clouds near active galactic nuclei, and planetary nebulae. We apply our models to five 
planetary nebulae and conclude that only BD + 30°3639 shows evidence of enhanced H 2 emission due to 
(high) nonequilibrium H 2 abundances. 

Subject headings: infrared: ISM: lines and bands — ISM: clouds — ISM: molecules 
planetary nebulae: general — ultraviolet: ISM 


1. INTRODUCTION 

The study of photodissociation regions (PDRs) is essential 
in understanding the properties of the interstellar medium. 
PDRs are neutral regions where far-ultraviolet (FUV) 
(6 eV < hv < 13.6 eV) photons (generally from OB stars or 
white dwarfs) play a significant role in the heating or chemistry 
of the gas. By this definition, PDRs include most of the mass of 
the interstellar medium, not only the atomic but also the 
molecular regions, where all the hydrogen is in H 2 and all the 
carbon is in CO but where residual FUV photons dominate 
the gas heating or the photodissociation of other important 
molecules, such as 0 2 or H 2 0. PDRs produce much of the 
continuum and line infrared emission in most normal or star- 
burst galaxies, for it is here that much of the starlight is 
absorbed by dust and where a fraction of the starlight energy is 
transformed to gas heating (Tielens & Hollenbach 1985, here- 
after TH85; Sternberg & Daigarno 1989; Hollenbach, Taka- 
hashi, & Tielens 1991, hereafter HTT91). Genzel (1991, 1992) 
recently reviewed the extensive observations of PDRs in the 
Milky Way as well as in many other galaxies. 

Most theoretical models of PDRs have focused on equi- 
librium chemical solutions (see Hollenbach & Tielens 1995 for 
a recent review of theoretical models). However, there has been 
a small number of time-dependent PDR calculations discussed 
in the last two decades. Hill & Hollenbach (1978) examined the 
propagation of the H/H 2 dissociation front when an O star 


suddenly turns on in a molecular cloud of density n = 10 3 -10 4 
cm ‘ 3 , where n = n H + 2n H2 is the hydrogen nucleus density. A 
prime concern was whether the shock associated with the 
expanding H ii region propagates into molecular or atomic 
gas. London (1978) studied the nonequilibrium H 2 front estab- 
lished when, for example, an OB star moves through a molecu- 
lar cloud. Roger & Dewdney (1992) extended the work of Hill 
& Hollenbach by including heating and cooling and expanding 
the parameter space to include densities of 10 ! -10 4 cm 3 and a 
range of stellar effective temperatures. Goldshmidt & Stern- 
berg (1995, hereafter GS95) recently studied the one- 
dimensional H/H 2 dissociation wave propagating through a 
semi-infinite molecular slab that is suddenly exposed to an 
FUV flux G 0 . Their focus was on the time-dependent fluores- 
cent H 2 emission produced for densities n= 10 2 -10 4 cm” 3 
and G 0 = 10°-10 6 , in units of the local interstellar FUV flux. 

Time-dependent PDR calculations are interesting because 
the timescales to establish H/H 2 equilibrium, r cq — 5 x 10 8 
n 1 yr, are often long compared with the timescales for varia- 
tion in the FUV field or the gas density. For example, we are 
primarily motivated by the PDRs created in the dense (n < 10 6 
cm 3 ) neutral shells around planetary nebulae. Here the 
central star may “turn on” its FUV luminosity in - 10 2 -10 3 
yr, and the density and FUV flux change on the shell dynami- 
cal timescale of -10 3 yr. Other applications include star- 
forming regions and clouds near active galactic nuclei (AGNs), 
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where clumpy material can be suddenly exposed to strong 
FUV radiation as clumps move out of shadowed regions. 

We present here a parameter study of constant-density, 
n = 10 3 -10 6 cm' 3 , semi-infinite slabs of molecular gas 
exposed at t = 0 to an FUV flux G 0 = 10 3 -10 6 . The focus is on 
the time-dependent photodissociation and fluorescence of H 2 
and the heating that accompanies the rapid FUV pumping and 
dissociation of the H 2 molecules. The relatively high densities 
lead to appreciable collisional de-excitation of the FUV- 
pumped, vibrationally excited H 2 , and the heating leads, in 
several cases, to a substantial thermal component to the H 2 
1-0 S(l) and 2-1 S( 1) intensities, thereby modifying the pure- 
fluorescent spectrum. We also calculate the atomic hydrogen 
column densities and the intensities of other observable lines 
that show interesting time dependence, including rotational 
transitions of OH and H 2 0 and the forbidden [Fe n] 1.64 /im, 
[O i] 6300 A, and [S ii] 6730 A. This work complements the 
results of GS95 by extending the parameter space to higher 
densities and by including a solution to the thermal balance 
and the line emission from a number of species. 

The paper is organized as follows: In § 2 we describe the 
physical processes incorporated into the time-dependent PDR 
code. In § 3 we present detailed results for a specific case 
(n = 10 s cm" 3 and G 0 = 10 5 ) and discuss simple analytic 
approximations to the evolution of the PDR. We also present 
and discuss the atomic hydrogen column densities and the 
emergent line intensities as a function of n, G 0 , and t. We briefly 
discuss possible applications of our general models in § 4 and 
summarize our conclusions in § 5. We discuss the numerical 
methods in the Appendix. We have incorporated this time- 
dependent PDR code in a much more comprehensive code 
designed to follow the evolution of expanding neutral shells or 
clumps surrounding planetary nebulae; these results will be 
presented in a subsequent paper (Natta & Hollenbach 1995). 

2. PHYSICAL PROCESSES AND NUMERICAL METHODS 
2.1. Chemistry and Timescales 

The neutral hydrogen is initially entirely molecular, but the 
onset of the FUV flux sends a wave of dissociation into the 
neutral layer. The time-dependent equation for the molecular 
hydrogen density n Hl is 

~ = */««„ - G 0 l 0 f{N H2 )e — , (1) 

where R f — 3 x 10~ 17 cm 3 s~ l is the rate coefficient for H 2 
formation on grains, n H and n Hl are the atomic and molecular 
hydrogen densities, / 0 = 4 x 10“ 1 1 s“ 1 is the unattenuated H 2 
photodissociation rate in the local interstellar radiation field, 
G 0 is the FUV photon flux incident on the slab in units of the 
local interstellar radiation field (about 2 x 10 7 photons cm’ 2 
s" 1 for the range 11 eV < hv < 13.6 eV), N and N H2 are the 
hydrogen nucleus and molecular hydrogen column densities, 
f(N H2 ) is the self-shielding factor for H 2 (GS95), and <r FUV = 1.4 
x 10“ 21 cm" 2 is the effective FUV grain-absorption cross 
section per hydrogen nucleus (TH85; GS95). In equation (1), 
the first term is the formation rate of H 2 on grains per unit 
volume; the second term is the FUV photodissociation rate of 
H 2 per unit volume. At sufficiently high molecular column 
densities, N H2 ;> 10 15 cm' 2 , FUV absorption in the Lorentz- 
ian wings of the Lyman-Werner lines initiates dissociation. 
The self-shielding factor then becomes f(N H2 ) = 
where p ~ 10 6 cm" 1 (Jura 1974; J. Black 1995, private 


communication). This condition generally applies for most 
regions of interest in the dissociation wave, as we shall discuss 
below. The characteristic timescale for H 2 to reach equilibrium 
is given, from equation (1), by (see also GS95) 


t 


eq 


1 

2 R f n 


~ 500 yr 



( 2 ) 


It should be noted that r eq is the timescale to reach equilibrium 
deep in the slab, where the fractional abundances of H and H 2 
are comparable (x H ~ x Hz ) and the FUV photodissociation is 
heavily shielded by dust and/or hydrogen molecules. As we 
shall discuss below, if this occurs at N > ffpuv ~~ 10 21 cm 2 , 
then this is also the timescale for the H 2 fluorescent intensities 
to drop from high, nonequilibrium values to much lower, equi- 
librium values. However, in the relatively unshielded regions 
where the H 2 equilibrium abundance x 2 q 1, the equilibrium 
timescale is much shorter than f eq and scales inversely with G 0 . 
We emphasize that our focus is on the parameter space 
G 0 /n > 10” 2 cm 3 , so that the equilibrium abundance of H 2 is 
very low until N > 10 21 cm" 2 , when dust shielding aids self- 
shielding and results in primarily molecular gas. 

In order to streamline the PDR code so that it may be used 
in our more complicated planetary nebula code (Natta & 
Hollenbach 1995), we have made simple approximations for 
the chemistry of other potential gas coolants. Our main goal is 
to calculate the emergent intensities of the important cooling 
transitions C ii 158 fx m, O i 63 /im, Si n 35 /mi, Fe ii 1.28, 1.64 
/an, O i 6300 A, and S ii 6730 A and the OH and H 2 0 rotation- 
al transitions, all of which may be important in the relatively 
warm and unshielded region N < 4 x 10 21 cm' 2 , where the 
H 2 emission originates. We have assumed gas-phase abun- 
dances of the elements x c ~ 3 x 10 ” 4 , x G = 5 x 10 4 , x Si = 
3.5 x 10“ 5 ,x s = 10' 5 , and x Fe = 10“ 6 . With the exceptions of 
OH and H 2 0, we have used the equilibrium PDR results of 
TH85 and HTT91 to estimate the abundances x(C + ), x(CO), 
x(O), x(Si + ), x(Fe + ), and x(S + ) as a function of column density 
N. We justify below the use of equilibrium abundances for 
these species, describe our simple fit to the equilibrium results, 
and discuss the calculation of OH and H 2 0 abundances. 

The timescale for CO to photodissociate and for carbon to 
ionize to C + in the neutral PDR is given by the photo- 
dissociation time of CO, which is prolonged by the self- 
shielding of CO but is nevertheless much shorter than the 
timescale for H 2 dissociation because JV(CO) < N H2 . For 
hydrogen column densities N > 10 21 cm" 2 , dust shielding sets 
in and CO persists. The maximum timescale required for 
C/C + /CO equilibrium is therefore given by the self-shielded 
photodissociation rate, calculated at N ~ 10 21 cm" 2 with 
N H2 ^ 5 x 10 20 cm" 2 and N co ~ 3 x 10 17 cm" 2 (all carbon 
in CO and all hydrogen in H 2 from the surface to N = 10 21 
cm' 2 ). Using these maximal-shielding column densities (the 
H 2 also partially shields the CO), we derive, using shielded 
rates from van Dishoeck & Black (1988), the carbon timescale 


*.**>*(%). 

This timescale is about 2 orders of magnitude [roughly 
(N C o/Nh 2 ) 1/2 ] smaller than the equivalent timescale for H 2 
photodissociation at the same column density. The timescale 
for other oxygen-bearing molecules, such as 0 2 , OH, and 
H 2 0, to photodissociate is even shorter since self-shielding is 
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unimportant: 

-„=:K>-yr(f). (4) 

These timescales are sufficiently short compared to the H 2 
timescales that we assume equilibrium PDR solutions apply 
for CO, C + ,and O. 

The equilibrium PDR results of TH85 and HTT91 show 
that for G 0 > 10 2 , n > 10 4 cm -3 , and GJn > 10 3 cm 3 , the 
carbon is ionized to hydrogen column densities of x 10 21 
cm -2 and the oxygen is entirely atomic to this same column 
density. Beyond this column density the carbon is in CO, and 
all oxygen not in CO remains predominantly atomic to N > 
10 22 cm -2 . Similarly, sulfur is ionized to IV ^ 4 x 10 21 cm -2 , 
while iron is ionized to greater depths. We therefore assume 
this simple prescription for the abundances of these species. 

However, inside the column density N < 4 x 10 21 cm -2 , a 
small but significant fraction of the oxygen may be converted 
to OH and H 2 0 by neutral-neutral reactions that proceed 
rapidly in hot molecular gas such as is produced when the 
dissociation wave propagates through the low column density 
gas. Although we find that the resultant enhanced abundances 
of OH and H 2 0 do not dominate the cooling in the cases 
considered, we include the emission from these species as pos- 
sible diagnostics of the time-dependent PDR. The important 
reactions include 

O + H 2 - OH + H , OH + H 2 -► H 2 0 + H . 

The rate coefficients for these reactions and the reverse reac- 
tions are given, for example, in Hollenbach & McKee (1989). 
We include in our PDR code this reaction sequence and the 
FUV photodissociation of OH and H 2 0 (see Hollenbach & 
McKee 1989 for FUV cross sections). 

This approximate chemistry underestimates the production 
of OH and H 2 0 because it ignores the ion-molecule reactions 
that can also produce enhanced OH and H 2 0 abundances in 
regions where the H 2 is higher than equilibrium. In addition, 
the enhanced OH abundance leads to an enhanced (though 
still small) abundance of CO. The CO chemistry is sufficiently 
complicated that we do not calculate the (warm) CO abun- 
dances in the enhanced-OH region, and therefore we cannot 
predict the emergent CO spectrum in this paper. However, the 
ion-molecule enhancements to OH and H 2 0 are small com- 
pared to their production by neutral-neutral reactions in the 
warm ( T > 300 K) regions, and the FUV photodissociation 
prevents CO from being the dominant coolant in these regions, 
so the overall temperature structure of the warm gas is well 
approximated. Our goal here is to calculate the H 2 and atomic 
(C It 158 pm, O i 63 ^m, Fe ii 1.2, 1.6 /*m, O i 6300 A, and S II 
6730 A) intensities and to make (lower limit) estimates of the 
OH and H 2 0 emission. 

2.2. Thermal Balance and Radiative Transfer 

The neutral gas is heated primarily by grain photoelectric 
heating (Bakes & Tielens 1994) and the heating caused by FUV 
pumping and formation pumping of the hydrogen molecule, 
followed by collisional de-excitation of the excited vibrational 
levels (see TH85; Hollenbach & McKee 1989). H 2 photo- 
dissociation also is a heating agent (TH85). These processes are 
well known and are often applied in the context of equilibrium 
chemistry; the novel aspect of their application to time- 
dependent PDRs is that, for example, nonequilibrium FUV 


pumping may dominate the heating as the dissociation wave 
sweeps through the molecular gas. 

Significant coolants in the neutral region include vibrational 
and rotational transitions of H 2 , collisional dissociation of H 2 , 
rotational transitions of CO, OH, and H 2 0, fine-structure 
transitions of O and C + , and grain cooling of the gas. We also 
include collisional excitation of Lya, [O i] 6300 A, [S n] 6730 
A, and [Fe n] 1.28, 1.64 jan, which are appreciable when the 
neutral region is driven to T > 3 x 10 3 K. The cooling rates 
for these processes have been taken from Hollenbach & 
McKee (1979, 1989) and TH85. The O i 63 /an, C ii 158 /an, 
CO, OH, and H 2 0 lines are treated with the escape- 
probability formalism described in Hollenbach & McKee 
(1979) because their line opacities can be large and self- 
absorption followed by collisional de-excitation may be impor- 
tant. A cooling time can be estimated by assuming LTE 
cooling by H 2 : 

~ 4Tj 7/2 yr , (5) 

where T 3 = T/1000 K and we have fitted the numerical cooling 
results of Hollenbach & McKee (1979) with a power law valid 
over the range 0.3 < T 3 <; 3. This timescale is sufficiently short 
that we assume thermal equilibrium while we calculate the 
nonequilibrium evolution of the hydrogen chemistry. On the 
other hand, the dissociation front moves rapidly through the 
cloud and causes a brief temperature spike in the gas at the 
dissociation front. We have checked the timescales and find 
that thermal equilibrium remains a reasonable assumption, 
even in the spike, for most cases of interest. In any event, there 
is little contribution to the overall observable line intensities 
from the hot spike, as discussed in § 3.3. 

The emergent intensities in the lines are found by integrating 
the escaping line emissivity through the thickness of the slab 
and dividing by 47i. For the case of the vibrational emission 
from H 2 , we add to the thermal component a “ nonthermaP 
contribution to the 1-0 5(1) and 2-1 5(1) intensities caused by 
the FUV pumping and formation pumping of the vibrational 
states. At high densities, n > n cr , the pumping contributions 
are decreased by a factor n cr /n, where n cr is the characteristic 
critical density for the collisional de-excitation of an excited 
vibrational level (see Burton, Hollenbach, & Tielens 1990 for 
details of this procedure). 

3. MODEL RESULTS 

3.1. The Time Evolution of the Dissociation Front 

Figures 1 and 2 present some detailed results for the evolu- 
tion of a slab with G 0 = 10 5 and n — 10 5 cm 3 . We use this 
case to demonstrate the evolution of the H 2 dissociation front 
and to examine in detail the temperature structure and the 
dominant heating and cooling processes in the gas. 

The evolution of x H2 as a function of N and t is shown in the 
bottom panel of Figure 1. Note that V, the hydrogen nucleus 
column density, is linearly related to depth l into the slab, 

/ = N/n. At early times, a steep dissociation front is apparent, 
where x H2 rises from its low equilibrium value ( < 1 0 3 ) to high, 
substantially molecular values (^>0.1) at a column density N*. 
We define N* to be the column density where the gas is half- 
atomic (x H = 0.5) and half-molecular (x H2 = 0.25). For these 
values of G 0 and n, N* < 10 21 cm 2 for t < 100 yr and grain 
attenuation is negligible at the front. Figure 1 shows that for 
these early times N* oc f 2 , or the velocity of the front v D oc 
dN*/dt oc t. For t > 100 yr, the front reaches N* > oyjy ~ 
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Fig. 1. — Temperature T {top) and H 2 fractional abundance x H2 ( bottom ) 
plotted as functions of the gas column density in the slab at different times. The 
slab has G 0 = 10 5 , n = 10 5 cm -3 . The depth into the slab / is also shown on 
the top scale. The steep dissociation front initially accelerates (v D oc f) into the 
slab. Enhanced column densities of warm H 2 are evident at early times. 


Log 1 (cm) 
15 16 



Fig. 2.— Heating and cooling rates of the dominant processes shown as 
functions of the gas column density at t = 100 yr for a slab with G 0 = 10 5 , 
n = 10 5 cm’ 3 . The upper panel plots heating rates; the solid line is the heating 
due to FUV pumping of H 2 followed by coilisional de-excitation, the dashed 
line the heating due to photoelectric effect on grains. The tower panel plots 
cooling rates; the solid line shows the rate of cooling due to H 2 rotational and 
vibrational line emission, the short-dashed line the rate of cooling due to 
coilisional dissociation of H 2 , the long-dashed line the O i 63 /im cooling, and 
the dotted line the cooling due to CO lines. 


10 21 cm 2 , grain attenuation lowers the FUV flux, and the 
front slows and eventually stops as equilibrium is reached at 
t > 10 3 yr. 

The steepness of the front is caused by “runaway” disso- 
ciation at small column densities. The self-shielding causes the 
highest dissociation rates to occur at the smallest N. As mol- 
ecules are dissociated at low N, the self-shielding decreases 
further, so the (already high) dissociation rates appreciably 
increase with time at low N , whereas at high N the shielding 
column densities have not substantially changed (N Hl ~ 0.5N) 
and the dissociation rates are constant in time. 

Ignoring dust extinction (i.e., N < 10 21 cm -2 ), one can ana- 
lytically solve equation (1) and derive the N and t dependence 
of x H and x„ 2 for the region N > N* since, at these columns, 
N H2 ^0.5N: 

x H = 2 ll2 fiI 0 G 0 N il2 t , (6) 

x„ 2 = 0.5(1 - x H ) , (7) 

valid for 

N* < N < o Fuv ~ 10 21 cm -2 . 

N* < N is the equivalent to the condition 

t KV^pioGoN^y 1 ~ ISON lit Gq^ yr, 

where N 21 = AT/10 21 cm -2 and G 05 = G o /10 5 . Equations (6) 
and (7) also assume that N > 10 15 cm 2 , so the self-shielding 
factor /(N„ 2 ) = p(0.5N)~ 1,2 , where ft = 10 6 cm -1 (see § 2.1). 
Either equation (6) or equation (7) can be used to derive N*{t), 
since the front tracks the position x H2 ~ 0.25 or x H ~ 0.5: 

N*(t) ~ 8/? 2 /q Gq t 2 ^ 1.2 x 10 21 Gq5 t\ cm 2 , (8) 

where t 2 = t/ 10 2 yr. We note that N* is greater than the 
column density where the front becomes extremely steep (see 
Fig. 1); this is because N* is defined as the column density to 
x H ~ 0.5, whereas the front steepens at x H ~ 0.1. The front 
velocity is calculated from equation (8) as 

v D ~ 16/? 2 /o G q 1 f — 77 Gqs Hs l h s ~ 1 , (9) 

where n 5 = n/10 5 cm 3 . Equation (8) and (9) are valid for 
N* < (TpJv - 10 21 cm 2 or, equivalently, t < t a , where t a is the 
time for the front to reach N — crpjy : 

t a =(8^ 2 / 2 G 2 a FUV ) 1/2 ^ 79G05 1 yr . (10) 

Equations (6)-(9) also implicitly assume that x H2 is signifi- 
cantly greater than its equilibrium value in the region of 
interest. This is always true for N < 10 21 cm -2 and t < t a if 
G 0 /n > 10“ 2 cm 3 (see§ 2.1). 

3.2. Comparison of v D with Previous Results 

Equation (9) seems to contradict the result in GS 95, where 
the “dissociation front” had constant speed with time. 
However, the discrepancy is only caused by a difference in the 
definition of the dissociation front. In GS95 the column density 
N to the front is identified as N the total (i.e., integrated to 
infinite depth) column density of atomic hydrogen at time t. By 
our definition, the column density to the front is N*, the 
column density at which the local abundances are half- 
molecular and half-atomic. These two column densities are 
different and have different time dependencies. Note from 
equation (6) that N H ocx H A 1/2 , so Afg 1 is dominated by the 
highest column density gas, where dust does not appreciably 
attenuate the FUV (in this case N ^ ^fuv ~ 10 21 2 f° r 
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larger column densities grain attenuation of the FUV rapidly 
lowers x H ). In other words, if we use Figure 1 at t = 10 yr as an 
example, it is not the completely dissociated N* ~ 10 19 cm -2 
that dominates VJJ 1 but the almost completely molecular gas at 
N ~ 10 21 cm -2 , where x H ~0.1 and N H — 7 x 10 19 cm 2 . 
Therefore, we are in agreement with GS95 that V„ l oc t at early 
times. From equation (6), taking the maximum N to be <rfj v , 
we obtain 

Njj 1 ~ 2 312 pi o Go ff Fuv 2 f • (It) 

With the GS95 definition of the dissociation front, the front 
velocity is proportional to dN^/dt, which is constant in time. 
However, for reasons made clear by Figure 1, we prefer our 
definition of the dissociation front, and we obtain v D x t for 
early times, t < t a (or N* < 10 21 cm -2 ), and for x H2 at N* 
significantly greater than the equilibrium value. 

3.3. Thermal Balance in Time-dependent PDRs 

Figure 2 illustrates some typical features of the thermal 
balance for the case n = 10 5 cm 3 and G 0 = 10 5 . The rapid 
rate per unit volume of FUV pumping of H 2 at the front, where 
x H2 is much larger than equilibrium and where there is little 
self-shielding, causes a temperature spike, as is evidenced by 
Figure 1. Figure 2 shows that, for the case studied, the heating 
in the spike is mainly caused by FUV pumping of H 2 followed 
by collisional de-excitation. The cooling is by collisional disso- 
ciation of H 2 and collisional excitation of ro-vibrational infra- 
red lines of H 2 . 

In the completely photodissociated region N < V*, grain 
photoelectric heating and H 2 FUV pump-heating are balanced 
by O i 63 pm cooling to produce a temperature of about 1500 
K in the atomic gas. Figure 1 also shows that the temperature 
at N ~ 10 21 cm -2 rises with time, which is caused by the 
decrease in the abundance of the coolant H 2 . At high column 
densities, IV>4 x 10 21 cm -2 , the grain photoelectric heating 
suddenly drops (Fig. 2) because the approximate carbon chem- 
istry produces a sudden drop in the electron abundance in the 
gas (which raises the positive grain charge and makes photoe- 
jection of electrons more difficult). Cooling is dominated here 
by O r 63 /an and rotational transitions of CO. 

3.4. FUV Pumping of H 2 , H 2 Fluorescence , and 
Thermal H 2 Emission 

GS95 extensively discuss H 2 fluorescence in low-density 
{n < n cr ~ 10 4 cm -3 ), time-dependent PDRs, and we briefly 
paraphrase their results here, with a somewhat different 
approach to explaining the limiting solutions. We first describe 
the low-density results, n < n cr , where collisional de-excitation 
of vibrationally excited H 2 is negligible and a “ pure- 
fluorescent ” spectrum emerges. 

Consider fluorescence when N Hl ~ 0.5V at N = ~ 

10 21 cm -2 . This case applies at all times for GJn < 10" 2 cm 3 
since self-shielding of H 2 enables the equilibrium abundance of 
H 2 to be x„ 2 ~ 0.5 at N — 10 21 cm -2 when the FUV-to- 
density ratio is this low. It also applies at early times, t < t a (or 
N* < 10 21 cm -2 ), for time-evolving PDRs with GJn > 10 - 2 
cm 3 . The column pump rate P FUV ( or the pump rate per unit 
PDR surface area) is about 9 times the column photo- 
dissociation rate of H 2 , or about 4.5 times the production rate 
of H : 

dN Xot 

P FUV — 4.5 — 4.5(2 3/2 /?/ 0 G 0 0>Jv 2 ) , (12) 


where we have made use of equation (11). The total intensity 
/, R in the infrared lines in this case is given by 

I\k = h = f>FU '' ~ 0.4G O5 ergs cm -2 s“‘ sr 1 , (13) 

47C 

where AE p is the average ro-vibrational energy delivered per 
FUV pump (~2.6 eV). Equation (13) can be easily understood. 
Essentially, H 2 wins the competition with dust for FUV 
photons, and most of the photons in the energy range 1 1-13.6 
eV are absorbed by H 2 and give rise to excitation and fluores- 
cence. In pure fluorescence the intensity of the 1-0 S(l) line is 
about 2% of /, R (see, e.g., GS95). The main result is that /, R is 
constant ( = /J at all times for G 0 /n < 10 -2 cm 3 and at early 
times, r < for GJn > 10" 2 cm 3 . Like the column density of 
atomic hydrogen, most of the line intensity originates from 
N ~ crpJv ~ 10 21 cm -2 . 

Now consider later times for the case GJn > 10 2 cm 3 . The 
gas becomes atomic at N ~ 10 21 cm -2 . For t a < t < f eq , the 
dissociation front slows as N* exceeds apjv and the grains 
attenuate the FUV. Most of the dissociation and, therefore, 
most of the pumps occur between the front N* and N* + <7 FU V 
Because of the exponential attenuation of dust, N* ~ N {J* for 
N > a Fl ! v . It follows that 

dN toi dN* f 

/ . R ° C_ df aC ”dr 0C J ! oGoe* FVvN f(N H2 )n H2 dl . (14) 

The sharp dissociation front makes the term n H2 e aFVyfN peak 
sharply between N* and N* + <t f J v so that 


dN* 

— oc I 0 G 0 e-™*\ (15) 

with the solution 

/|R ° c ^ _oc ^ oGot ' 1 

for the case GJn £ 10 - 2 cm 3 and t a <t<t c ^. 

Finally, we consider the pumping rate in equilibrium, t > 
t cq . GS95 show that in this case J IR oc N 5J* ~ trfjy, ignoring a 
logarithmic term of order unity. This can be understood by 
noting that in equilibrium there are nine FUV pumps per H 2 
formation and that the total column formation rate of H 2 is 
R/ hVh*. Thus, we obtain 


/lR = /«, ^ 


9R r nNu 



To summarize these results for n < n CT : if GJn < 10“ 2 cm 3 , 
then / IR = l ff , with l a given by equation (13). If GJn > 10 2 
cm 3 , then /, R = I a for t < t 9 , / IR ~ l a {tjt) for t„ < t < r eq , and 
JiR - ^ eq t > t eq , with t eq given by equation (2), t a by equa- 
tion (10), and / eq by equation (17). 

Figure 3 illustrates these low-density results for GJn > 10 2 
cm 3 . The cumulative H 2 intensity (the intensity produced by 
gas between 0 and N) is plotted as a function of N at various 
times. The low-density case n — 10 3 cm -3 is plotted so that the 
results are not complicated by the effects of collisional de- 
excitation. For t < t a ~ 3 x 10 3 yr, note that most of the inten- 
sity is produced by gas at N ~ 10 21 cm -2 and that the total 
integrated intensity is constant in time. At later times, t a < t < 
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Fig. 3.— Cumulative fluorescent 1-0 5(1) H 2 intensity (i.e., the intensity 
produced by gas between 0 and N) shown as a function of the gas column 
density N at various times for a case where G 0 = 10 4 , n = 10 3 cm 3 . The figure 
shows that at early times most of the intensity is produced at N ^ <7pJ v ^ 10 21 
cm -2 . 


f eq ~ 5 x 10 5 yr, the intensity is produced somewhat deeper 
into the slab at the (attenuated) dissociation front, and / 1R x 
t~ l . Finally, in equlibrium, the intensity is mainly produced in 
the final logarithmic interval of AN, where H rapidly converts 
to H 2 (around N - 3 x 10 21 cm“ 2 ). 

At high density, n > n cn the FUV fluorescent component of 
the H 2 emission is reduced by a factor n/n CT . This energy goes 
into heat, and Figures 1 and 2 show that that this can raise the 
temperature of considerable column densities of H 2 to 
T > 1000 K. In this case H 2 can act both as the dominant heat 
source and as the dominant coolant. The heat is delivered by 
the collisional de-excitation of FUV-pumped vibrational 
levels. The cooling arises from the collisional excitation, and 
subsequent radiative decay, of low-lying vibrational levels. 
Therefore, much of the fluorescent pump energy emerges as 
thermal emission from low-lying H 2 vibrational transitions. To 
add to the complication, n Cf is temperature sensitive and is 
smaller at higher temperatures, thereby enhancing collisional 
de-excitation. 

3.5. Model Results of Line Intensities 

Figures 4, 5, and 6 show the main results of significance to 
observers. They plot the total (thermal plus fluorescent) emer- 
gent intensity in the H 2 1-0 5(1) line, the fluorescent com- 
ponent of the H 2 1-0 5(1) line, and the ratio of the total 2-1 
S(l)/l-0 5(1) line intensities. Figure 4 demonstrates the depen- 
dences discussed above. For low values of G 0 /n, less than 10 2 
cm 3 , there is no time dependence of / fR ; / 1R — / CT . On the other 
hand, for low values of n, n < n cr , and high values of G 0 /h, 
greater than 10 2 cm 3 , we see the following sequence for I lR {t): 
first / ff , then and finally / eq . The main result shown in 

Figure 4 is that the H 2 intensities can show strong, time- 
dependent variations to times as long as 5 x 10 8 /n yr. The 


nonequilibrium intensities can be at least an order of magni- 
tude larger than the equilibrium values. 

Figure 5 shows the effects of collisional de-excitation of the 
pump-excited vibrational levels. At low density, n = 10 3 cm 3 , 
2-1 5(l)/l-0 5(1) — 0.5, the pure-fluorescent value (with a 
slight variation at t = 100 yr for G 0 > 10 5 due to the tem- 
perature dependence of n CT ). There is also a pure-fluorescent 
ratio for densities as high as 10 6 cm 3 if G 0 < 10 3 since the 
pump region is cold and the critical densities are high (> 10 6 
cm -3 ) for such cold temperatures. However, for G 0 > 10 4 and 
n> 10 3 -10 4 cm' 3 , the 2-1 S(l)/l-0 5( 1) ratio rapidly 
approaches a thermal value of order 0.1. This ratio has often 
been interpreted to indicate shock emission, but Figure 5 
demonstrates that collisionally altered FUV pumping can 
produce this ratio. Figure 5 also shows that the timescale for 
the 2-1 S(l)/l-0 5(1) ratio to reach a steady state is somewhat 
less than t cq . 

Figure 4 and 6 also show the effects of collisional de- 
excitation of FUV-pumped H 2 . Figure 6 plots the fluorescent 
component of the 1-0 5(1) intensity in order to provide an 
estimate of the intensities of the higher vibrational lines (i> > 2). 
The higher vibrational states will not have a significant 
thermal component because they lie A E/k > 18,000 K above 
the ground state and cannot be significantly excited by colli- 
sions with gas at T ~ 1000 K. Therefore, their intensities can 
be estimated from Figure 6 by using the pure-fluorescent ratios 
of the lines in question to 1-0 5(1) (see, e.g., Sternberg 1988) 
and multiplying the ratio by the 1-0 S(l) intensity in Figure 6. 

Figure 7 plots the time-dependent intensities in the strong 
(surface brightness greater than 1 0 “ 6 ergs cm ~ 2 s " 1 sr 1 ) cool- 
ants OH, H 2 0, Fe ii 1.64 \i m, O i 6300 A, and S u 6730 A, 
which show some interesting dependences. Individual lines 
from these coolants have surface brightnesses greater than 
10“ 6 ergscm“ 2 s“ 1 sr -1 only for the case n = 10 6 cm“ 3 ,so we 
have only plotted this case. Because only the high-density cases 
have detectable lines, steady state is rapidly reached in 10-100 
yr. The strong lines of C II 158 /mi, O i 63 /mi, and Si n 35 /an 
are not shown because they show little time dependence and 
steady state values can be taken from TH85 and HTT91. 
Therefore, Figure 7 demonstrates that the emission from other 
species will probably not help diagnose nonequilibrium in 
time-dependent PDRs. Observers can take the values from 
equilibrium PDR models to compare with observed line inten- 
sities. 

3.6. Model Results of the Evolution of the Atomic 
Hydrogen Column 

The H i 21 cm line can be used to trace the total atomic 
hydrogen column N^ x in planetary nebulae and other poten- 
tially time-dependent PDRs. Likkel et al. (1992) present upper 
limits to the 21 cm fluxes from two compact planetary nebulae 
(BD +30°3639 and AFGL 618) and review previous 21 cm 
studies of planetary nebulae. Taylor, Gussie, & Pottasch (1990) 
have observed H i 21 cm in a few planetary nebulae and con- 
clude that the H i mass detected cannot adequately account for 
the mass expected to have been lost by the asymptotic giant 
branch precursor star. It is possible the mass is in H 2 and that 
the dissociation front has not yet propagated through the bulk 
of the mass. 

Figure 8 therefore plots N[f as a function of u, G 0 , and t. The 
time-dependent creation of atomic hydrogen is clearly seen, 
and the final equilibrium columns are generally of order 10 21 
cm -2 because dust attenuation of the FUV preserves the 
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Fig. 4. — Total intensity of the H 2 1-0 S{1) line as a function of time; the units are ergs cm " 2 s“ 1 sr“ \ Each panel refers to a different value ofG 0 , as labeled. The 
solid lines are for n = 10 3 cm -3 , the dotted lines for n = !0 4 cm 3 , the short-dashed lines for n = 10 5 cm“\ and the long-dashed lines for n = 10 6 cm 3 . In the 
G 0 = 10 6 panel, log n is explicitly labeled on the lines. For GJn > 10" 2 cm 3 , the intensity drops with time from a peak value I a to / eq (see text). 


hydrogen molecules at higher columns. For low values of G 0 /n , 
<10“ 2 cm 3 , the equilibrium values of are significantly 
smaller than 10 21 cm -2 because of H 2 self-shielding. 

4. APPLICATIONS 

GS95 discuss applications of time-dependent PDR models 
to reflection nebulae and to neutral gas in starburst galaxy 
nuclei or in active galactic nuclei. They also discuss application 
to planetary nebulae (PNs), and we amplify that discussion 
here since PNs are our primary motivation in developing this 
code. The neutral densities in PNs often exceed 10 4 cm 3 
(otherwise the gas is rapidly ionized), and the timescales for 
changes in the FUV flux and gas density are short, of order 
100-1000 yr. The PDR region lies in a thin shell of thickness 
A R outside of the ionized nebula at a distance R from the 
central white dwarf. Because A R/R <£ 1, the one-dimensional 
slab models well approximate the conditions in the PDRs in 
PNs (Natta & Hollenbach 1995). 

Huggins (1992, 1993) reviews the status of observations of 
neutral gas associated with PNs. Over 33 PNs have been 
detected in H 2 2 /im vibrational emission. Although pure- 
fluorescent emission has been observed (for example, in 


Hubble 12 and BD + 30°3639; Dinerstein et al. 1988), in many 
cases the 2-1 5(l)/l-0 5(1) ratio is of order — 0.1 (e.g., Zucker- 
man & Gatley 1988). This has often been interpreted as arising 
from shock emission, perhaps caused by the ejection of a dense 
shell that supersonically overtakes the slowly expanding red 
giant wind of the PN progenitor. However, as has been already 
argued (Hollenbach 1988; Sternberg & Dalgarno 1989; Burton 
et al. 1990), FUV pumping at high density can produce similar 
“thermal ” ratios. 

Without detailed models of individual PNs (which is the 
subject of Natta & Hollenbach 1995), it is difficult to argue that 
the observed H 2 emission may reflect FUV pumping of non- 
equilibrium H 2 , even if the incident FUV flux can be estimated. 
As Figure 4 demonstrates, for fixed G 0 , a given intensity can be 
produced at early times in low-density, nonequilibrium PDRs 
or at late times from equilibrium PDRs of higher density. GS95 
argue that the densities in the planetaries NGC 7027 and BD 
+ 30 3639 are sufficiently low(<3 x 10 4 cm 3 ) that nonequi- 
librium models are needed to explain the observed intensities. 
However, in the case of NGC 7027, the 2-1 5(l)/l-0 5(1) ratio 
is low, —0.07, and not the -0.5 expected from the pure fluo- 
rescence of a low-density PDR (Tanaka et al. 1989). In addi- 
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Fig. 5— Ratio of the 2-1 S(l) to the 1-0 S(I) line intensities as a function of time. Lines denote values of n as in Fig. 4. A ratio of 0.5 indicates pure fluorescence; 
lower ratios indicate contributions from thermal collisions. 


tion, a PDR model fit to the O i 63 /tm and C n 158 fim 
intensities predicts PDR densities of 7 x 10 4 cm" 3 (Ellis & 
Werner 1985) and an FUV flux of G 0 = 4 x 10 4 (GS95). The 
peak 1 — 0 5(1) intensity is about 1 x 10" 3 ergs cm -2 s" 1 sr" 1 
(Graham et al. 1993a, b; Cox et al. 1995). Figures 4 and 5 show 
that a reasonable fit to the 1-0 5(1) and the 2-1 5(1) line 
intensities can be obtained if G 0 ~ 10 5 and n ~ 10 5 cm" 3 , con- 
sistent with the fine-structure lines. However, the fit requires 
t > 10 3 yr and nearly equilibrium conditions. Therefore, NGC 
7027 may not show evidence of nonequilibrium chemistry. 

BD +30°3639 is a better case because the spectrum looks 
like pure fluorescence, which suggests low densities, n < 10 4 
cm" 3 , yet the 1-0 5(1) peak intensity is also about 1 x 10" 3 
ergs cm" 2 s" 1 sr" 1 (Graham et al. 1993a). GS95 quote G 0 ~ 2 
x 10 4 . Figures 4 and 5 show that a fluorescent 2-1 5(1)/ 1-0 
5(1) ratio of 0.5 and the observed 1-0 5(1) intensity are produc- 
ed if n < 10 4 cm" 3 and if t < 10 3 yr. The nonequilibrium H 2 
intensity in this case is at least of an order of magnitude higher 
than the equilibrium value. 

Zuckerman & Gatley (1988) observed H 2 emission in three 
planetaries, NGC 2346, NGC 6720 (the Ring), and NGC 6853 
(the Dumbbell). The peak 1-0 5(1) intensities are about 


1 x 10" 4 ergs cm 2 s" 1 sr" 1 and the inferred G 0 < 10 3 in all 
three cases. The 2-1 S(l)/l-0 S(l) ratios are all less than 0.1. 
The neutral densities are poorly constrained, but an estimate 
can be made with the assumption that the neutral gas is in 
pressure equilibrium with the ionized gas. Zuckerman & 
Gatley estimate electron densities in the ionized gas of 300- 
1000 cm" 3 , which indicate neutral densities of n ~ 10 4 cm" 3 . 
They argue that shock emission, with H 2 excited by collisions 
with other hydrogen molecules, will be too weak to explain the 
observed intensities, so a PDR model is worth considering. 
Figure 4 shows that the observed peak intensities can (barely) 
be obtained in nonequilibrium models, if one assumes the most 
optimistic (high) values of G 0 . If the neutral densities are < 10 4 
cm" 3 , then nonequilibrium models apply and t < 10 4 yr is 
required. However, Figure 5 predicts that all models with 
G 0 ^ 10 3 and n < 10 6 cm" 3 have fluorescent ratios 2-1 5(1)/ 
1-0 5(1) ~ 0.5, contrary to the observations. We conclude that 
even a time-dependent PDR model has difficulty explaining 
the observations. Perhaps the best solution at this time is one 
proposed by Zuckerman & Gatley: the H 2 is excited by colli- 
sions with hydrogen atoms (which have much larger rate coef- 
ficients than hydrogen molecules) in shock waves. 
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Fig. 6.— Same as Fig. 4, but for the fluorescent component of the 1-0 S(l) line. This intensity may be multiplied by the pure-fluorescent ratios of high vibrational 
(t?) lines to 1-0 S( 1 ) as in, e.g., Sternberg ( 1988) to obtain predictions on the time-dependent intensities from high-r lines that have no thermal component. 


5. SUMMARY AND CONCLUSION 

In this paper the time dependence of photodissociation 
regions was examined specifically for the case of a fully molecu- 
lar slab that is suddenly exposed to an FUV flux G 0 . Our 
time-dependent PDR code follows the photodissociation and 
fluorescence of H 2 and includes the calculation of the heating 
that accompanies the rapid FUV pumping and dissociation of 
the hydrogen molecule. Both thermal collisions and fluores- 
cence contribute to the excitation of H 2 vibrational lines. If 
G 0 > 10 4 and n > 10 5 cm -3 , thermal collisions dominate the 

1- 0 S(l) and 2-1 S(l) line intensities at all times; thermal colli- 
sions contribute significantly even for densities as low as 
n ~ 10 4 cm -3 if t < 10 3 yr. These thermal collisions produce 

2- 1 S(l)/l-0 S(!) ratios <0.1. Such ratios have often been 
taken to indicate shock excitation; however, in this case it is 
simply FUV-heated gas. Figures 4, 5, and 6 present a param- 
eter study in n, G 0 , and t of the 1-0 S{\) and 2-1 S(l) line 
intensities and a method to calculate the emergent intensities of 
higher vibrational (pure fluorescent) transitions. 

As the dissociation wave sweeps into the slab, the H 2 vibra- 
tional intensities are time dependent only if G 0 /n > 10 2 cm 3 . 
In this case, the time dependence proceeds as follows: initially, 


for t < t a (see eq. [10]), the intensity has a constant value of 
(eq. [13]); for t 9 < t < f eq (eq. [2]), the intensity drops as 
finally, for t > f eq , the intensity asymptotically 
approaches the constant / (eq. [17]). The dissociation wave is 
characterized by a steep “front,” which initially accelerates 
through the slab with velocity proportional to time. The front 
slows once it penetrates to a column density N ~~ cr F J v ~ 10 21 
cm -2 , where dust attenuation begins to sharply reduce the 
dissociating-photon flux. At all times, most of the H 2 
vibrational-line intensities are produced at N ~~ 10 21 cm 2 , 
and it is the drop in H 2 abundance at this column density that 
causes the H 2 intensity to drop. Our calculations extend to 
higher density the recent results of GS95 for the fluorescent H 2 
emission. 

We also presented (Figure 7) the surface brightnesses of O i 
6300 A, S ii 6730 A, Fe n 1.64 /im, and the total rotational 
contributions of OH and H 2 0. The intensities of these tran- 
sitions also show interesting time dependence and are observ- 
able if n > 10 6 cm -3 . The atomic lines also require G 0 > 10 5 
for observability. The total atomic hydrogen column densities 
are shown in Figure 8 for possible application to H i 21 cm 
observations. 
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Fig. 7. — Intensity of some metal and molecular lines as a function of time. Long-dashed lines refer to models with G 0 — 10 3 , short -dashed lines to G 0 = 10 4 , 
dotted lines to G 0 = 10 5 } and solid lines to G 0 = 10 6 . In all cases, n = 10 6 cm “ 3 . The OH and H 2 0 panels refer to total intensities in all pure-rotational lines of v = 0. 
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APPENDIX 

NUMERICAL METHOD 

We assume a constant-density, semi-infinite slab, initially fully molecular and then illuminated at time t = 0 by an FUV flux G 0 . 
Although the temperature varies in the slab, we assume that turbulence dominates the pressure so that the density is always 
constant and uniform (as did the original equilibrium PDR models of TH85). The code numerically integrates equation (1), using 
logarithmic intervals of hydrogen nucleus column density and time. The initial column density bin has a column density AN — 10 17 
cm 2 , and the initial time step is At = 1 yr. Thereafter, we use typically 200 column density bins per decade and 100 timesteps per 
decade. We have reduced the bin sizes and the time steps by factors of 2 until successive solutions match to better than a few percent. 
We have also compared the results with analytic solutions for the time dependence (see § 3) and the equilibrium solution. Finally, we 
have run several cases with constant temperature and low density to compare the results with the values quoted by GS95. We find 
good agreement (differences of less than ±30%, and even less when we adjust for the slightly different values of/ 0 and<r FUV in our 
code). 

The only difficulty in this integration is the large difference between the timescale for dissociation in the unshielded FUV field, 
t min ^ (10 3 /G o ) yr, compared to the timescale to reach final equilibrium, r eq ^ (5 x 10 8 /n) yr. If we are forced to take time steps of 
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Fig. 8— Total atomic hydrogen column density N™ plotted as a function of t, G 0 , and n. The different lines denote different densities ; lines denote values of n as in 
Fig. 4. 


less than f min , then we would be forced to take 10 8 time steps (on each of 1000 spatial grid points) in order to follow the evolution of 
the case n = 10 3 cm 3 and G 0 = 10 5 , for example! We therefore use the following algorithm, which allows time steps much greater 
than f min . We integrate equation (7) over a small time step At and find, at a column density N, 

x Hl (t + At) = [x„ 2 (t) - x eq > A,/,ch + x eq , ( A 1 ) 


where x„ 2 is the H 2 abundance and where 


= R f nt ch 


(A2) 


is the equilibrium abundance of H 2 , 


f ch = [2g / n + G 0 /o/(N2 2 k' , '™ vN ] 1 


(A3) 


is the characteristic timescale for change of the H 2 abundance, and 

= [N„ 2 (t) + N H2 (t + At)]/2 <A4) 

is the average shielding column density between t and t + At. The column density N„ 2 (f + At) is known because we start at the 
surface and solve for x„ 2 (f + At) from the surface into column density N. With this method, even if At P t ch f mj „, the solution 
moves in an appropriate way toward the equilibrium solution. 
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